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Abstract 

We examine the properties of a recently proposed model for antigenic variation in malaria 
"which incorporates multiple epitopes and both long-lasting and transient immune responses. We 
show that in the case of a vanishing decay rate for the long-lasting immune response, the system 
exhibits the so-called "bifurcations without parameters" due to the existence of a hypersurface 
of equilibria in the phase space. When the decay rate of the long-lasting immune response is 
different from zero, the hypersurface of equilibria degenerates, and a multitude of other steady 
states are born, many of which are related by a permutation symmetry of the system. The 
robustness of the fully symmetric state of the system was investigated by means of numerical 
computation of transverse Lyapunov exponents. The results of this exercise indicate that for a 
vanishing decay of long-lasting immune response, the fully symmetric state is not robust in the 
substantial part of the parameter space, and instead all variants develop their own temporal 
dynamics contributing to the overall time evolution. At the same time, if the decay rate of 
the long-lasting immune response is increased, the fully symmetric state can become robust 
provided the growth rate of the long-lasting immune response is rapid. 

1 Introduction 



Several pathogens, including Plasmodium falciparum malaria and African trypanosomes, achieve 
immune escape by the so-called antigenic variation (see a recent review by Gupta [7]). The lat- 
ter essentially refers to a process by which a pathogen keeps changing its surface proteins, thus 
preventing antibodies from recognizing and destroying it. Antigenic variation is achieved by ex- 
ploiting a large repertoire of antigenic variants that differ in some of their epitopes. An important 
requirement here is that the variants must not be expressed all at the same time, as otherwise the 
resulting immune response will detect and destroy all of them, thus terminating the infection. 

"This work was partially supported by the ATRJVVO grant from the James Martin 21st Century School, Uni- 
versity of Oxford. 
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Here we examine a particular model of antigenic variation for P. falciparum malaria, put forward 
by Recker et al. [ID]. Within this framework, each variant is assigned one major epitope, which is 
unique to that variant, and also several minor epitopes that are shared between different variants. 
Both types of epitopes elicit epitope-specific responses, but in the case of the minor epitopes these 
are cross-protective between variants that share them. A critical feature of the model is that the 
immune response to the major epitope (uniquely variant-specific) is long-lasting in comparison with 
the immune responses (frequently cross-protective) to the minor epitopes. Under these conditions, 
the dynamics may be characterized by sequential domination of different variants. Thus, the 
conclusion of the model is that effectively the host immune system can itself be responsible for 
prolonging the malaria infection and causing chronicity. Through numerical simulations and by 
analysing a caricature of the model involving complete synchrony between variants, Recker and 
Gupta [11] have shown that stronger cross-protective immune responses lead to prolonged length 
of infection and reduced severity of the disease. This was explained by the conflicting interaction 
of cross-protective and variant-specific immune responses. 

In this paper we perform a detailed study of this model with particular emphasis on stability 
aspects, as well as possible bifurcation scenarios. First, we consider the case when the variant- 
specific immune responses to the major epitopes do not decay. In this case, the phase space of 
the system possesses a very peculiar geometry with a high-dimensional surface of equilibria having 
different types of stability. We show, it is exactly this curious structure that causes successive re- 
appearance of different malaria variants in the dynamics until the specific immune responses reach 
sufficient protective level to prevent further appearance of given variants in the dynamics. If the 
specific immune responses can decay (even slightly), the dynamics is qualitatively different, as the 
phase space geometry changes significantly. Now it contains a large number of distinct equilibria 
with different number of non-zero variants, some of which can be related by the permutation 
symmetry of the system. 

An interesting tool for investigating the dynamics is the imposition of synchrony among the 
variants. From mathematical perspective, in the case of complete synchrony the dimension of the 
system is drastically reduced. Here, we use the tools of synchronization theory to investigate the 
robustness of such state. 

The outline of this paper is as follows. In the next section the model of cross-reactive immune 
response to malaria is introduced and its basic properties are discussed. Section 3 contains the 
analysis of a particular case when the decay rate of a long-lasting immune response vanishes. 
Numerical simulations will be presented that illustrate the behaviour of the system in this case. 
A general situation of arbitrary non-decaying specific immune responses is considered in Section 
4. In section 5 the stability of the fully symmetric state of the system is investigated by means of 
numerical computation of transverse Lyapunov exponents. The paper concludes in section 6 with 
a discussion. 

2 Model definition 

In this section we use the above-mentioned multiple epitope description to introduce a model of 
the interaction of malaria variants with the host immune system. Our derivation follows that of 
Recker et al. [TDJ with some refinements. 

It is assumed that each antigenic variant i consists of a single unique major epitope, that elicits 
a long-lived (specific) immune response, and also of several minor epitopes that are not unique to 
the variant. Assuming that all variants have the same net growth rate eft, their temporal dynamics 
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is described by the equation 

yi((j) - azi - a'wi), (1) 



dt 

where a and a' denote the rates of variant destruction by the long-lasting immune response z% and 
by the transient immune response Wi, respectively, and index i spans all possible variants. The 
dynamics of the variant-specific immune response can be written in its simplest form as 

dz- 

= fyi - fizi, (2) 

with f3 being the proliferation rate and \i being the decay rate of the immune response. Finally, 
the transient (cross-reactive) immune response can be described by the minor modification of the 
above equation (J2j) : 

~df =P'J2 y i -m'w*. (3) 

where the sum is taken over all variants sharing the epitopes with the variant yi. We shall use 
the terms long-lasting and specific immune response interchangeably, likewise for transient and 
cross-reactive. 

To formalize the above construction, one can introduce the adjacency matrix A, whose entries 
Aij are equal to one if the variants i and j share some of their minor epitopes and equal to zero 
otherwise. Obviously, the matrix A is always a symmetric matrix. Prior to constructing this matrix 
it is important to introduce a certain ordering of the variants according to their epitopes. For this 
purpose we shall use the lexicographic ordering, as explained below. To illustrate this, suppose 
we have a system of two minor epitopes with two variants in each epitope, which is the simplest 
non-trivial system of epitope variants. In this case, the total number of variants is four, and they 
are enumerated as 



1 11 

2 12 

3 21 

4 22 



(4) 



It is clear that for a system of m minor epitopes with n\ variants in each epitope, the total number 
of variants is given by 

m 

N = Hm. (5) 

1=1 

Now that the ordering of variants has been fixed, it is an easy exercise to construct the adjacency 
matrix A of variant interactions. For the particular system of variants (j3J), this matrix has the form 



(6) 



In general, for a system of two minor epitopes with m variants in the first epitope and n variants 
in the second, the matrix A will be an mn x mn block matrix consisting ofnxn blocks of ones 
along the main diagonal, with the rest of the matrix being filled with n x n identity matrices. 
For simplicity, in the rest of the paper we will concentrate on the case of two minor epitopes, but 
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Figure 1: Interaction of malaria variants in the case of two minor epitopes with two variants in 
each epitope. 



the results can easily be generalized for arbitrary number of minor epitopes. Using the adjacency 
matrix one can rewrite the system (pQ)-([3]) in a vector form 




y(4>l N ~ az 
P'Ay - /x'w, 



a'w) 



(7) 



where y = (1/1,2/2, ■■■,Un) etc., Ijv denotes a vector of the length ./V with all components equal to 
one, and in the right-hand side of the first equation multiplication is taken to be entry-wise so that 
the output is a vector again. 

To better understand the symmetry of the system it is convenient to represent graphically 
relations between different variants. Figured] shows such relations in the case of two minor epitopes. 
One can observe that within each horizontal and each vertical stratum, the network of variants is 
characterized by an "all-to-all" coupling [3]. Besides this, if the number of variants in both minor 
epitopes is the same, then there is an additional reflectional symmetry. Formally this means the 
system is equivariant with respect to the following symmetry group [6] 



r> J S m x S n , 
5= S x a 



m 7^ n, 



S m x S m x Z2, m = n. 



(8) 



This construction can be generalized in a straightforward way for a larger number of minor epitopes. 
It is noteworthy that each of the variants has exactly the same number of connections to other 
variants. 

We finish this section by noting that the system ([7]) is well-posed, in that provided the initial 
conditions for this system are non-negative y(0) > 0,z(0) > 0,w(0) > 0, the solutions satisfy 
y(t) > 0,z(i) > 0,w(i) > 0, for all t > 0. 

Remark. In many cases it is reasonable to assume the initial conditions for the system ([7]) to be 
of the form (y > 0, z = 0, w = 0). A possible exception is when the immune system has already 
built-up a long-lasting response from prior exposure to a certain variant. In this case, the initial 
condition for the system ([Tj will contain non-zero entries for some of z variables 
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3 The case of non-decaying specific immune response 



We begin our analysis of the system ([7]) by considering a particular case of vanishing decay rate of 
the long-lasting immune response \x = (some partial results for this case have been obtained in 
|10j). In this case the only steady states of this system are given by 

y = w = 0, z = const. (9) 

This is a rather degenerate situation as the fixed points are not separated in the phase space, but 
rather form an TV-dimensional hypersurface with each point of it being a fixed point of the system 
([7]). Linearization near one such fixed point has the eigenvalues (— //) of multiplicity N, zero of 
multiplicity N, and the rest of the spectrum is given by 

(j) — azi, cj) — aZ2, <j) — aZtf. (10) 

The generalized eigenvectors of the zero eigenvalue correspond to the N directions along the hyper- 
surface of the fixed points. As long as there is at least one Zi < (fr/a, the corresponding steady state 
is a saddle, otherwise it is a stable node. From the dynamical systems perspective, the case \x = 
corresponds to the so-called bifurcations without parameters [5]. Indeed, in the space {y = w = 0} 
as one crosses the hyperplane Zj = 4>/a, one of the eigenvalues crosses zero along the real axis. 
Furthermore, since the hypersurface {y = w = 0} is in general high-dimensional, the cases of two or 
more eigenvalues crossing zero at the same time (this happens along the lines Z{ = Zj = cp/a,i ^ j) 
are still generic, and these lead to "bifurcations" of a higher co-dimension. It is important to 
note that all these bifurcations are of the steady state type and there is no possibility of a Hopf 
bifurcation that could lead to temporally periodic solutions. 

Figure [2] shows numerical simulations of a typical behaviour in the system ([7]) for [i = 0. 
These results were obtained by integrating the system (J7J) using the variable order solver based 
on backward differentiation formulas to account for the stiffness of the system [13]. Initially most 
variants have quite high amplitudes, but as the time progresses, their amplitudes decrease as 
illustrated in Fig. [2]^ a). Figures (b) and (d) illustrate this feature in more detail by showing 
the dynamics of a single variant and its specific immune response. With each subsequent re- 
appearance of the variant, the specific immune response to it is building up, and ultimately it 
reaches a protective level, which prevents this variant from ever re-appearing in the dynamics. As 
suggested by Fig. Efc), sometimes more than one variant appear at the same time, and this is very 
good from the immune system perspective, as it allows simultaneous destruction of all of these 
variants. The question of synchronization between different variants will be investigated in Section 
5. 

The fact that the system exhibits the jumps from one variant to another can be explained 
by the existence of the above-mentioned hypersurface of equilibria. When one of the variants 
decays, the trajectory approaches the neighbourhood of the hypersurface of equilibria, and since 
all points on this hypersurface are saddles of different dimensions, the trajectory is pushed away 
along the unstable manifold of one of these fixed points. This behaviour is reminiscent of that in 
the neighbourhood of a heteroclinic cycle [HE], with the major difference being that in the present 
case the nodes of the cycle are not distinct but rather form a smooth hypersurface. There is a 
clear separation of time scales in the dynamics: the trajectories move quickly to/away from the 
invariant w — z plane, and then they slowly move towards the hyper-axis y = w = before the next 
iteration. With time, the phase space excursions between subsequent returns to the equilibrium 
manifold become shorter (they are restricted by the ever growing z variables), and eventually all 
trajectories converge to a point T = (y = w = 0, z = ((ft/a)!^). Similar behaviour takes place in 
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Figure 2: Temporal dynamics of the system (|7|) with two variants in each of the two epitopes. 
Parameter values are a = a' = 10 _3 ,/3 = f3' = 10 _4 ,/i = 0, fj! = 0.02,0 = 7.5 [TU]. a) Time 
evolution of the variants yi. b) Dynamics of a single variant y\. c) The close-up of the initial stage 
of evolution of the variants, d) Dynamics of the long-lasting immune response z\. 



the phase coordinates of other variants, which all approach the point T. The importance of such a 
point for understanding the dynamics has been previously highlighted, for instance, in the analysis 
of adaptive control systems, where it gave rise to bad point bifurcations at which the close- loop 
systems could never be stabilised [12]. 

If during time evolution, the trajectory reaches the hyperplane z% = <p/a for some i, and at least 
one of yi or Wi is different from zero, then this trajectory will escape the basin of attraction of the 
point (y = w = 0, z = (^>/a)ljv) and instead it will asymptotically converge to the hypersurface 
of equilibria with the value of z% > (ft/ a without any further phase space excursions. This will 
happen provided the initial amount of a given variant is high enough. Figure [3] (a) shows in red 
the projection of the stable manifold of the point on the reduced phase space of a single variant 
together with a representative trajectory in blue. In the same figure a trajectory in green illustrates 
the scenario in which the protective level of immune response is reached within one parasitemia 
peak, and hence there are no further oscillations. In Fig. [3] (b) we show the close-up of the phase 
dynamics in the neighbourhood of the hypersurface of equilibria. One can clearly observe recurrent 
oscillations of parasitemia, during which the specific immune response is monotonically increasing 
until it reaches the protective level. 



6 



Next we would like to discuss the issues of peak dynamics and the threshold for chronicity, 
which have been previously studied in Recker and Gupta [11J. The chronicity threshold is defined 
as the critical ratio of the variant destruction rates 7 C = a' /a, such that if 7 < 7 C , then during 
the first peak the protective level of immunity will be reached, so that the system will display no 
further oscillations. There are several simplifying assumptions, which have to be made in order to 
derive analytical expressions for the solutions needed for the analysis of peak dynamics. First of 
all, it is reasonable to assume that all variants in the full system (|7|) are identical, and therefore 
this system can be replaced by 

y = y(4> — az — a'w), 

z = Py, 

w = P'n c y - fi'w, 

where n c is the number of connections for each variant, and y^ = y, etc. The second assumption is 
that for a single parasitemia peak the cross-reactive immune response does not have time to decay, 
i.e. for a peak dynamics we have // = 0. This reduces the system to 



y = y(4>- 
z = Py, 
w = (3'n c y 



az 



aw) 



(11) 



Assuming zero initial conditions 2(0) = w(0) = 0, which correspond to the absence of pre-existing 
specific or cross-reactive immune responses, the analytic expression for the solutions of the system 
(TTT1) can be found as 



«'«> = 7p 



+ \JlC\i) tanh 



(t + c 2 ) 



y(t) 



Ci 

P 



1 — tanh 



(t + c 2 



(12) 



where the integration constants C\ and C 2 are given by 



Ci = /3yo + 7T7, C 2 



arctanh 



'Ct-Pyo 



Ci 



and ip is defined as 



a/3 + a'nf3' 



Initially, y(t) monotonically increases, until it reaches its peak of y = C\ exactly at t = —C2, after 
which y(t) is monotonically decreasing. Due to the symmetry of the solution, at t = —2C 2 , y has 
the same value as it had at the initiation of parasitemia peak. By considering the equation for y(t), 
one can argue that if at the end of the parasitemia peak the combined specific and cross-reactive 
immune response has reached the protective level of (ft /a, then this will prevent further oscillations. 
Evaluating z and w at the end of parasitemia peak, we find the threshold for chronicity as 



z(-2C 2 ) + w(-2C 2 ) > 



la 



if 7 < 7c 



P 



a 

- = 2 + , 

a n c p' 



If 7 < 7c, then during the first peak y should be sufficiently high to allow the build-up of protective 
immunity. Conversely, if 7 > j c , then y will be too low for protective immunity to be reached 
within one peak, and therefore the system will display further oscillations |llj . 
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It is important to note that the trajectories shown in red and blue in Fig. [3] satisfy the condition 
for chronicity 7 < j c , but still for these trajectories the protective level of immune response is not 
reached within a single parasitemia peak. The reason for this discrepancy is due to the fact that 
in the system describing the peak dynamics, the cross-reactive immune response does not decay, 
because if it did, then at the end of the parasitemia peak the combined immunity would be below 
the protective level, and therefore further oscillations would occur, as shown in Fig. [3j This also 
highlights the importance of initial conditions, and in particular, the initial amount of variants, 
which may play a crucial role in whether or not the protective immunity level will be reached within 
one peak for the same parameter values. 



4 General case 

In the previous section we considered the case \i = 0, in which the long-lasting immune responses 
can only grow with time, unless they are saturated at the level preventing further re-emergence of 
particular variants. For /x > 0, the situation is drastically different as the hypersurface of equilibria 
no longer exists, and instead, it degenerates into two separate steady states. One of these is the 
origin (y = z = w = 0), which is always a saddle with an TV-dimensional unstable manifold and 
a 2iV-dimensional stable manifold. The other steady state originating from the hypersurface of 
equilibria is the fully symmetric equilibrium 

^ l aft fx' + a'n c j3' /i' 1 a/3{j,' + a'n c f3' fj,' 

(13) 

* a/3[i' + a'n c /3'fi ' 

where n c = m + n — 1 is the number of connections for each variant. Using Fig. 1, this number 
can easily be interpreted as the number of elements in the horizontal and vertical strata, to which 
the current variant belongs. When considered in the context of a reduced system (jlip . in which 
all variants are assumed to behave in the same manner (see next section for further analysis of 
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Figure 4: a) Boundary of Hopf bifurcation in a parameter space, b) The two-parameter Hopf 
boundaries for different values of the specific response decay rate \i = 0.005 (solid), \x = 0.01 
(dashed) and \x = 0.015 (dotted). In (a) the fully symmetric steady state is unstable above the 
coloured surface, and in (b) it is unstable above each of the corresponding curves, (c) Temporal 
dynamics of variants 1 (red), 2 (blue), 3 (green) and 4 (cyan), (d) Temporal dynamics of variant 1 
showing periodic oscillations with large intermittent period of quiescence. 



this case), this steady state is stable for all values of parameters, as shown in Recker and Gupta 
[TT] . At the same time, this result does not hold for the full system 0, as the stability of the 
fully symmetric equilibrium does depend on parameters of the system. More specifically, the fully 
symmetric equilibrium can undergo Hopf bifurcation, thus giving rise to periodic occurrences of 
parasitemia peaks. 

Figures H] (a) and (b) show the boundary of the Hopf bifurcation in the parameter space of the 
system ([7]) with two variants in each of the two minor epitopes. These figures indicate that the 
higher is the decay rate of variant specific immune response, the larger should be the values of the 
relative immune efficiency 7 = a' /a and that of a ratio of proliferation rates /3'/f3 to guarantee 
the occurrence of the Hopf bifurcation. The corresponding temporal evolution of variants in the 
parameter regime beyond the Hopf bifurcation is illustrated in figures (c) and (d). One can observe 
periodic oscillations of all four variants, which have approximately the same maximum amplitudes 
and are slightly out-of-phase with each other. The plot of the dynamics of one variant shown in 
Fig. H] (d) indicates that peaks of parasitemia corresponding to this variant have decreasing ampli- 
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tudes, and after several occurrences there are large periods of time when the variant is quiescent. 
Cross-reactivity between different variants causes subsequent re-appearance of large-amplitude os- 
cillations after such periods of quiescence. The reason for this is as follows. The long-lasting and 
cross-reactive immune responses show anti-phase oscillations, which are quite regular both in am- 
plitude and in period. These oscillations lead to a slightly irregular oscillations of the combined 
rate of variant destruction az + a'w. Intervals of parasitemia peaks correspond to the combined 
variant destruction rate oscillating around the critical value of (ft with long-lasting immune response 
increasing. After such intervals, the combined variant destruction rate stays above (ft keeping the 
variant absent from the dynamics, and during this time the long-lasting immune response wanes, 
until it starts to recover during the next cycle. We emphasize that this dynamics can only occur in 
the case when the long-lasting immune response can decay, hence this feature could not be observed 
in the previously analysed case of \i = 0. 

Besides the origin and a fully symmetric equilibrium, the system also possesses (2^ — 2) steady 
states characterized by a different number of non-zero variants yi. One should notice that the 
symmetry of the system mentioned earlier implies that for a given number of non-zero variants, 
many of the corresponding steady states are symmetry-related. At the same time, one can identify 
several clusters of the steady states with different values of the steady states which cannot be 
transformed into each other by a symmetry. For example, if we consider the system with two 
variants in each of the two minor epitopes, then there exist six steady states with two non-zero 
variants. Introducing the notation 

y = y = 

1 P'a'n + aPn n 2 2/3'a , /i + a/3/i'' 

the steady states with non-zero variants 12, 13, 24 and 34 form one cluster: 

E12 = (Y 2 ,Y 2 ,0,0,Z 2 ,Z 2 ,0,0,2W 2 ,2W 2 ,W 2 ,W 2 ), 
E 13 = (Y 2 ,0,Y 2 ,0,Z 2 ,0,Z 2 ,0,2W 2 ,W 2 ,2W 2 ,W 2 ), 
E 2i = (0, Y 2 ,0, Y 2 ,0, Z 2 , 0, Z 2 ,W 2 , 2W 2 , W 2 , 2W 2 ), 
E 24 = (0, 0, Y 2 , Y 2 , 0, 0, Z 2 , Z 2 , W 2 , W 2 , 2W 2 , 2W 2 ), 

while the steady states with non-zero variants 14 and 23 are in another cluster 

E u = (Y l5 0, 0, Yi, Z x , 0, 0, Z u Wx, 2W X , 2Wi, Wi), 
E 23 = (0, Y u Yi, 0, 0, Z u Z u 0, 2W U W u W u 2W X ), 

where Z\ j2 = f3Yi^/n and Wi )2 = f3'Yx >2 /fj,' . All the steady states in the first cluster are related 
by permutation, and the steady states in the second cluster are also related by some permutation, 
but the steady states from the first cluster cannot be related to those in the second cluster. The 
reason for this becomes clear if one more closely analyses the structure of the adjacency matrix A 
given in ([6|). In the case of a steady state E^ from the first cluster, both rows i and j of matrix 
A contain ones in positions i and j (i.e. the variants i and j cross-react with each other), while 
in the case of the second cluster the rows i and j contain only a single one in either position i or 
position j (i.e. the variants i and j are completely unrelated). Due to this difference the steady 
states from the two clusters are different and it is impossible to change from one cluster to another 
by permutation. 

As far as stability of the steady states different from the origin and the fully symmetric equi- 
librium is concerned, they all are saddles of different dimensions. Even though they are unstable 
as steady states, it is possible for some of them to form some sort of a heteroclinic cycle. 
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Remark. In the case when the number iV of malaria variants participating in the dynamics exceeds 
four, the symmetry of the system increases the co-dimension of the Hopf bifurcation for the fully 
symmetric steady state. Moreover, the purely imaginary eigenvalues at the Hopf bifurcation would 
coincide, thus creating extra complications for the analysis by virtue of increasing the dimension of 
the centre manifold. Some details of possible bifurcation scenarios in systems with an "all-to-ah" 
coupling can be found in [3j 0], and the extension of those results should provide an insight into 
the effects of symmetry on the dynamics of system ((7|). The complete analysis of these effects will 
be presented elsewhere. 



5 Robustness of the fully symmetric solution 

An interesting dynamical regime occurs when, by virtue of initial conditions or time evolution, the 
system behaves in such a way that all variants are indistinguishable from each other, in other words, 
the system is in a state of complete symmetry. In this case, the dimension of the system reduces 
drastically from 3iV to just three. As several insightful results have been obtained for this case |llj . 
it is important to study how robust this state of complete symmetry is with respect to perturbations 
that attempt to break the symmetry. To characterize stability properties of the symmetric state 
one can use transverse Lyapunov exponents, as is customary in the studies of synchronization, see, 
for instance [8]. By analogy with synchronization theory we shall call the hypersurface of complete 
symmetry a symmetry manifold. 

Writing {yi = y + y~i, z; L = z + Zi,Wi = w + Wi), one can split the total dynamics into that inside 
the symmetry manifold 

y = y((j) — az — a'w), 

z = (3y- nz, (14) 
w = (3'n c y - fi'w, 

and the linearized dynamics in the transverse direction given by 

y 

DF(yl N ,zl N ,wl N ) I z 




w 



4> — az — o.'w)In —otylN —ct'ylN 
DF(yl N ,zl N ,wl N ) = | pi N -yuljv On 



(15) 



Here n c is again the number of connections of a given variant, On and I/v denote N x N zero 
and unity matrices, respectively. The minimal condition for the stability (or robustness) of the 
symmetric state is that the maximum Lyapunov exponent associated with the system (|15p has 
negative real part [8]. By solving equations f)15|) in combination with ()14|) . we determine the 
dependence of the leading Lyapunov exponent on the system parameters. 

In Figure [5] we show the results of numerical simulations for the maximal transverse Lyapunov 
exponent A™ ax . In both plots we kept the rates of variant destruction equal to each other a = a' 
and also the proliferation rates were taken to be the same (5 = f} 1 . Figure [5] indicates that for small 
values of a = a' or f3 = (3', the fully symmetric state of the system is transversely unstable, as 
signified by the positive transverse Lyapunov exponent. This means that in such parameter regime 
different variants will not synchronize in time, hence it is unlikely to observe the fully symmetric 
state in experiment. However, when the variant destruction rates/proliferation rates are increased, 
the fully symmetric state becomes transversely stable, i.e. independently on initial conditions for 
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(a) (b) 




Figure 5: Maximal transverse Lyapunov exponent as a function of parameters, a) Parameter values 
are <fi = 7.5, /3 = /?' = 1, // = 0.02. b) Parameter values are <fi = 7.5, a = a' = 1, // = 0.02. 
Different colours represent different decay rates of the long-lasting immune response: /x = (red), 
fi = 0.01 (blue), fi = 0.02 (green) and fi = 0.03 (magenta). 

each particular variants, they will all ultimately follow the same time evolution. The increase in the 
decay rate of the specific immune response fi plays a stabilizing role, since it lowers the values of of 
the maximal transverse Lyapunov exponent. The robustness of the fully symmetric state appears 
to be independent on the relative efficiency 7 = a' /a of immune responses. 

6 Conclusions 

In this paper the temporal behaviour in a model of antigenic variation in malaria has been studied 
from a dynamical systems perspective. Using the model of immune response to multiple epitopes, 
we have demonstrated that when the long-lasting immune response does not decay, the system pos- 
sesses a high-dimensional surface of equilibria, and these exhibit steady-state bifurcation without 
parameters, i.e. some part of the surface of equilibria consists of saddles of different dimensions, 
while another part contains stable nodes. The existence of these two parts of the surface of equi- 
libria with different stability properties accounts for the observed patterns of behaviour of malaria 
variants, when different variants exhibit out-of-phase parasitemia peaks that decay with time. If 
the initial amounts of all variant are not very large, then phase space excursions between successive 
re- appearances of the variants become shorter as the time grows, and eventually all trajectories 
approach the single steady state T characterized by all coordinates equal to each other and equal 
to the value at the boundary between the saddles and the nodes on the surface of equilibria. If, 
however, an initial amount of a given variant is sufficiently high, then a trajectory with such ini- 
tial condition will escape the basin of attraction of the above-mentioned point T by reaching the 
protective level of long-lasting immune response to a given variant while having either a non-zero 
transient response to this variant or a non-zero amount of the variant itself. In this case, the 
eventual time evolution of the solution will be different in that it will also approach the surface of 
equilibria but now it will be above the critical protective level without any further excursions in 
the phase space. 

When both variant-specific and cross-reactive immune responses are allowed to decay with a 
certain rate, the dynamics are quite different. In this case the surface of equilibria disintegrates, 
and instead the phase space of the system contains a large number of distinct fixed points many 
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of which are related to each other by the permutation symmetry of the variants. At the same 
time, they may form separate clusters which are not related by symmetry. Provided the decay rate 
of the specific immune response is high enough, the fully symmetric equilibrium will exhibit Hopf 
bifurcation, thus giving rise to periodic oscillations of the variants. These oscillations appear to be 
out-of-phase for different variants, and such oscillations are separated by extended time intervals 
during which the amount of a variant is very small. 

In order to investigate to what extent the results obtained in the approximation of complete 
symmetry between variants describe the general patterns of behaviour, we have numerically com- 
puted the transverse Lyapunov exponents of the fully symmetric state. This analysis indicates that 
while the fully symmetric state is not robust to small perturbations for small proliferation/variant 
destruction rates, the robustness is restored as these rates increase. In this case the dynamics of 
the completely symmetric system faithfully represents that of the full original system. Finally, we 
note that the robustness of complete synchronization between variants increases with the decay 
rate of the specific immune response. 
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